{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e284e211",
   "metadata": {},
   "source": [
    "# Uncovering Associations among Genes and Phenotypes with SEAHORSE\n",
    "Author: Enakshi Saha<sup>1</sup>\n",
    "\n",
    "<sup>1</sup>Department of Biostatistics, Harvard T. H. Chan School of Public Health, Boston, MA, USA."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4d81a5e",
   "metadata": {},
   "source": [
    "# 1. Introduction\n",
    "\n",
    "SEAHORSE (Serendipity Engine Assaying Heterogenous Omics-Related Sampling Experiments) is a tool for exploratory data analysis designed to allow users to identify both expected and unexpected associations within large data sets; SEAHORSE allows users to ask and answer questions about individual variables or cohorts. Given gene expression and clinical data from a particular cohort such as a specific tissue from the Genotype Tissue Expression (GTEx) Project or a particular cancer type from The CancerGenome Atlas (TCGA), seahorse identifies pairs of most correlated genes, and also genes and biological pathways strongly associated to any given phenotype.\n",
    "\n",
    "You need to set the `runserver` parameter to 1, if you're running this vignette on the server. Otherwise, if the vignette is ran locally this parameter has to be set to 0.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26a663a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "runserver=1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18c9380e",
   "metadata": {},
   "source": [
    "If you are running this vignette locally, you need to install SEAHORSE through the netZooR package by following the instructions here: https://netzoo.github.io."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48eecd1c",
   "metadata": {},
   "source": [
    "Next, we need to load the `netZooR` package to use SEAHORSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "235e5829",
   "metadata": {},
   "outputs": [],
   "source": [
    "library(netZooR)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b22a8d",
   "metadata": {},
   "source": [
    "We also need to load the package fgsea to run gene set enrichment analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70446ab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "library(fgseaRes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c94e3df7",
   "metadata": {},
   "source": [
    "# 2. A Simple Example with Toy Data Generated from GTEx Lung Samples\n",
    "We download TPM-normalized RNAseq data from http://gtexportal.org/ and select 200 random lung tissue samples from the GTEx (version 8) database. We choose one categorical phenotype \"SEX\" and one numeric phenotype \"HGHT\" (height) to demonstrate how SEAHORSE identifies associations from gene expression and phenotypic data.\n",
    "\n",
    "The phenotype HGHT is a dbGap-protected and is not publicly available. To protect data privacy we add random noise to the HGHT variable to generate toy phenotypic data for each of the 200 individuals.\n",
    "\n",
    "We filter out lowly expressed genes (genes with count < 1TPM in at least 5% of all samples).\n",
    "\n",
    "The expression and the phenotypic data can be downloaded as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8288f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "if (runserver==0){\n",
    "    system(\"curl -O  https://netzoo.s3.us-east-2.amazonaws.com/netzoo/netZooR/tutorial_datasets/seahorse/GTEx_lung_expression_toydata.txt\")\n",
    "    system(\"curl -O  https://netzoo.s3.us-east-2.amazonaws.com/netzoo/netZooR/tutorial_datasets/seahorse/GTEx_lung_phenotype_toydata.txt\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7afb97a0",
   "metadata": {},
   "source": [
    "On the server, the expression and the phenotypic data can be loaded as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de95ba94",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = \"/opt/data/netZooR/seahorse/\"\n",
    "# Load expression data\n",
    "expression_data <- read.csv(paste0(file_path, \"GTEx_lung_expression_toydata.txt\"), sep=\"\")\n",
    "# Load phenotypic data\n",
    "phenotype_data <- read.csv(paste0(file_path, \"GTEx_lung_phenotype_toydata.txt\"), sep=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6caba70c",
   "metadata": {},
   "source": [
    "Let us first print the first 5 rows of the phenotypic data to see what kinds of phenotypes we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18fdccb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(phenotypic_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "031dbed1",
   "metadata": {},
   "source": [
    "We need to create a vector of same length as the number of phenotypes recording the type (\"numeric\" or \"categorical\") of each phenotypic variable. We see from above that \"SEX\" is a categorical variable with two categories and \"HGHT\" is a numeric variable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fee5c2d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "phenotype_dictionary = c(\"categorical\", \"numeric\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc040d8a",
   "metadata": {},
   "source": [
    "Next we load a list of biological pathways to perform gene set enrichment analysis for every phenotypic variable, using the association between each gene and the phenotype under consideration. Here we use the KEGG biological pathways downloaded from https://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp?collection=H that contains 187 pathways. The pathway file can be downloaded as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a29f71f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if (runserver==0){\n",
    "    system(\"curl -O  https://netzoo.s3.us-east-2.amazonaws.com/netzoo/netZooR/tutorial_datasets/seahorse/c2.cp.kegg.v2022.1.Hs.symbols.gmt\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c03e196",
   "metadata": {},
   "source": [
    "On the server, the list of pathways can be loaded as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0325de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load KEGG pathways\n",
    "pathways_KEGG <- gmtPathways(paste0(file_path, \"c2.cp.kegg.v2022.1.Hs.symbols.gmt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec569eac",
   "metadata": {},
   "source": [
    "Finally, we run SEAHORSE using the gene expression and the phenotypic data, the phenotype dictionary vector we created and the list of pathways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8406593",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = seahorse(expression = expression_data, phenotype = phenotype_data, \n",
    "                   phenotype_dictionary, pathways = pathways_KEGG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77209f20",
   "metadata": {},
   "source": [
    "The SEAHORSE output contains three components: \n",
    "- the coexpression matrix containing Pearson's correlation for every pair of genes.\n",
    "- the measure of association between every gene and every phenotype.\n",
    "- a table containing the results from gene set enrichment analysis (GSEA) for every phenotype.\n",
    "\n",
    "Let us examine each component one by one.\n",
    "\n",
    "First, the gene-gene coexpression matrix can be used to identify most correlated genes. For example, suppose we are interested to find top 10 genes most strongly correlated (positive or negative correlation) with the gene \"EGFR\". The mutated forms of the EGFR gene have been found in some types of cancer, including non-small cell lung cancer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab82cf0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Coexpression matrix\n",
    "cor = results$coexpression\n",
    "# Coexpression of all genes with EGFR\n",
    "cor_EGFR = cor[which(rownames(cor) == \"EGFR\"),]\n",
    "# Top 10 genes (excluding EGFR itself) most correlated with EGFR\n",
    "cor_EGFR[order(abs(cor_EGFR),decreasing = T)[2:11]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f39b0c62",
   "metadata": {},
   "source": [
    "We can identify top biological pathways that are most differentially expressed between males and females in lung tissue using the association component of the SEAHORSE results. The association between a categorical phenotype and gene expression is measured by the p-value of an ANOVA test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "982cafbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation between height (HGHT) and all genes\n",
    "cor_sex = results$phenotype_association[[\"SEX\"]]\n",
    "# Top 10 genes most differentially expressed by sex\n",
    "cor_sex[order(cor_sex)][1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac500b68",
   "metadata": {},
   "source": [
    "Next, we identify top genes that are most (positively or negatively) associated with height using the association component of the SEAHORSE results. The association between a numeric phenotype and gene expression is measured by the Pearson correlation coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75a09a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation between height (HGHT) and all genes\n",
    "cor_hght = results$phenotype_association[[\"HGHT\"]]\n",
    "# Top 10 genes most correlated with height (HGHT)\n",
    "cor_hght[order(abs(cor_hght), decreasing = T)][1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "419d3ac3",
   "metadata": {},
   "source": [
    "We can also identify the biological pathways significantly correlated (say, at a false discovery rate cutoff 0.05) with height in the lung tissue. The following code prints top 5 biological pathways most associated with height."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "654e39ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene set enrichment analysis results for height (HGHT)\n",
    "GSEA_hght = results$GSEA[[\"HGHT\"]]\n",
    "# Top 5 biological pathways most associated with height.\n",
    "GSEA_hght[order(GSEA_hght$padj),][1:5,]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b405808",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "1- Ben Guebila, Marouen, Tian Wang, Camila M. Lopes-Ramos, Viola Fanfani, Des Weighill, Rebekka Burkholz, Daniel Schlauch et al. \"The Network Zoo: a multilingual package for the inference and analysis of gene regulatory networks.\" Genome Biology 24, no. 1 (2023): 45.\n",
    "\n",
    "2- Lonsdale, John, Jeffrey Thomas, Mike Salvatore, Rebecca Phillips, Edmund Lo, Saboor Shad, Richard Hasz et al. \"The genotype-tissue expression (GTEx) project.\" Nature genetics 45, no. 6 (2013): 580-585.\n",
    "\n",
    "3- Data Coordinating Center Burton Robert 67 Jensen Mark A 53 Kahn Ari 53 Pihl Todd 53 Pot David 53 Wan Yunhu 53, and Tissue Source Site Levine Douglas A 68. \"The cancer genome atlas pan-cancer analysis project.\" Nature genetics 45, no. 10 (2013): 1113-1120.\n",
    "\n",
    "4- Kanehisa, Minoru. \"The KEGG database.\" In ‘In Silico’Simulation of Biological Processes: Novartis Foundation Symposium 247, vol. 247, pp. 91-103. Chichester, UK: John Wiley & Sons, Ltd, 2002.\n",
    "\n",
    "5- Mailman, Matthew D., Michael Feolo, Yumi Jin, Masato Kimura, Kimberly Tryka, Rinat Bagoutdinov, Luning Hao et al. \"The NCBI dbGaP database of genotypes and phenotypes.\" Nature genetics 39, no. 10 (2007): 1181-1186."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
